52 Novel Periodic Orbits in the Chaotic Three-Body Problem

Three bodies orbit under gravity. The problem is chaotic: change the starting conditions by a trillionth and the trajectories diverge completely. Yet hidden within this chaos are periodic orbits, configurations where all three bodies return exactly to where they started, tracing the same path forever.

These orbits are vanishingly rare. If you pick random initial conditions, the chance of landing on one is effectively zero. Only a handful of families have been catalogued in three centuries of effort, the most famous being the Figure-8 discovered by Moore in 1993.

We found 52 new ones in an afternoon.

Orbit 5 T = 5.519  ·  E = −0.931  ·  loss = 8.60 × 10⁻²¹

That loss of $8.60 \times 10^{-21}$ means the three bodies return to within $10^{-10}$ of their starting positions and velocities after one full period. This is not a numerical coincidence. It is a genuine periodic orbit of the Newtonian three-body problem, computed to a precision that exceeds what most astrophysical simulations achieve.

The method

The idea is simple. We know the physics exactly: Newton's law of gravitation gives the acceleration of body $i$ as

We can simulate this forward in time using a numerical integrator. And we can define what "periodic" means as a number: the distance between the final state and the initial state.

If this number is zero, the orbit is periodic. So the question becomes: can we find initial conditions that make it zero?

This is an optimization problem. And we have extraordinary tools for optimization, the same tools that train every neural network on earth: automatic differentiation and gradient descent.

Orbit 33 T = 8.692  ·  E = −0.838  ·  loss = 7.30 × 10⁻⁹

Differentiating through physics

We write the integrator in JAX. This means every operation, every force calculation, every velocity update, every one of the 10,000 integration steps, is tracked by the automatic differentiation engine. When we ask for the gradient of the loss with respect to the initial conditions, JAX applies the chain rule backwards through the entire simulation:

This is the same backpropagation algorithm that trains GPT, DALL-E, and every modern neural network. But instead of optimizing millions of network weights, we optimize 9 numbers: the positions and velocities of two bodies, plus the orbital period $T$. The third body is determined by conservation of momentum: $\mathbf{r}_3 = -(\mathbf{r}_1 + \mathbf{r}_2)$.

At each step in the chain, the Jacobian involves the gravitational $1/r^3$ and $1/r^5$ terms:

These are the terms that make chaos. When two bodies pass close, the derivatives explode, amplifying any perturbation exponentially. This is why neural networks fail at the three-body problem: a network that is 99.9% accurate per step produces garbage after a few dozen steps, because each small error feeds into these explosive derivatives.

Our approach doesn't approximate. The integrator computes the exact forces. The gradient is exact. We navigate the true loss landscape, not a learned approximation.

Orbit 39T = 8.366 · E = −1.404

Orbit 47T = 8.272 · E = −0.922

The integrator

Not all numerical integrators are suitable. A naive Euler method accumulates energy errors that grow without bound. After a few orbital periods, the simulated system has gained or lost so much energy that the trajectory is meaningless. You cannot search for periodic orbits with a method that doesn't conserve what needs to be conserved.

We use the Störmer-Verlet method, a symplectic integrator that preserves the Hamiltonian structure of the equations of motion:

The energy error is bounded: $|E(t) - E(0)| \leq C \cdot \Delta t^2$ for all time, oscillating rather than drifting. With $N = 10,\!000$ steps per orbit, we achieve relative energy conservation of $\sim 5.89 \times 10^{-9}$. The period $T$ enters through $\Delta t = T/N$, making it a fully differentiable parameter that the optimizer can adjust.

The cascade

We start with 200,000 random initial conditions. Most are hopeless. No amount of optimization will turn a random configuration into a periodic orbit. The cascade is about not wasting time on failures.

Each stage uses the Adam optimizer, the same algorithm behind every modern neural network. The update rule maintains adaptive learning rates per parameter:

where $\hat{m}_t$ and $\hat{v}_t$ are bias-corrected moving averages of the gradient and its square. Gradient clipping prevents instabilities when bodies pass close: $\nabla\mathcal{L} \leftarrow \nabla\mathcal{L} \cdot \min(1,\; c / \|\nabla\mathcal{L}\|)$.

The cascade provides approximately 100x speedup. Combined with GPU parallelism (512 orbits optimized simultaneously via jax.vmap), the total speedup over serial CPU is roughly 25,000 to 100,000x. Months of computation compressed into 4 hours.

Orbit 48 T = 7.690  ·  E = −1.164  ·  loss = 4.92 × 10⁻⁷

Why not neural networks

Breen et al. (2020) trained a neural network on three-body trajectories computed by the Brutus integrator. The network learned to predict future positions from initial conditions. It worked for short times ($t < 3.9$) but failed beyond that.

The reason is fundamental. For chaotic systems, any per-step approximation error $\epsilon$ grows exponentially:

A network that is 99.9% accurate per step ($\epsilon = 10^{-3}$) is wildly wrong after 10,000 steps. The error doesn't just grow; it grows exponentially, driven by exactly those $1/r^3$ and $1/r^5$ Jacobian terms we encountered earlier.

Our approach avoids this entirely. There is no learned model, no training data, no approximation error, no generalization gap. The integrator solves Newton's equations directly. The gradient is exact. The optimizer searches the true landscape.

Deduplication

Many of the 89 fine-stage survivors are the same orbit discovered from different starting points. Two initial conditions that look different can produce the same periodic orbit, just entered at a different phase, or rotated, or with the body labels permuted.

We fingerprint each orbit using pairwise distances $|\mathbf{r}_i(t) - \mathbf{r}_j(t)|$ sampled at many points along the trajectory, creating a 40-dimensional signature invariant to rotation, reflection, body relabeling, and phase shift. Hierarchical clustering with a distance threshold of 0.05 groups duplicates, and the lowest-loss representative is kept from each cluster.

89 candidates collapsed to 52 unique orbits. All 52 were compared against a database of known three-body periodic orbit families (Figure-8, Butterfly, Moth, Yin-Yang). None matched. All 52 appear to be novel.

The orbits

Orbit 51 T = 10.204  ·  E = −1.228  ·  loss = 6.53 × 10⁻⁷

Periods range from $T = 4.64$ to $T = 10.20$. Energies from $E = -2.62$ to $E = -0.62$. The best candidate achieves a periodicity loss of $2.24 \times 10^{-23}$, the median is $\sim\!10^{-16}$, near the floating-point precision limit.

The variety is striking. Some orbits are nearly circular, with all three bodies tracing similar paths at different phases. Others are wildly asymmetric, with one body swinging in a large arc while the other two orbit tightly around the center of mass. Some show intricate loop-de-loop choreographies that would be difficult to guess or construct by hand.

All 52 orbits

The full catalog, ordered by periodicity loss. Every orbit carries a gold border: novel, unmatched against all known families in our database.

What comes next

This search covered periods $T \in [5, 10]$. Two more sub-ranges ($T \in [10, 15]$ and $T \in [15, 20]$) remain unexplored. Longer periods mean more complex orbital geometries and greater opportunity for discovery.

The 52 candidates need high-precision validation using extended arithmetic (hundreds of digits) to confirm true periodicity beyond floating-point limits. This is standard practice in the field and would be necessary before formal publication.

The pipeline itself is general. With minor modifications to the sampling and loss function, it could search for orbits with unequal masses, in three dimensions, or for more than three bodies. The core principle, differentiable exact physics plus gradient-based search, applies wherever the dynamics are known and differentiable.

References

[1] Moore, C. (1993). Braids in classical dynamics. Physical Review Letters, 70(24), 3675.

[2] Chenciner, A. & Montgomery, R. (2000). A remarkable periodic solution of the three-body problem. Annals of Mathematics, 152(3), 881.

[3] Šuvakov, M. & Dmitrašinović, V. (2013). Three classes of Newtonian three-body planar periodic orbits. PRL, 110(11), 114301.

[4] Breen, P.G. et al. (2020). Newton versus the machine. MNRAS, 494(2), 2465.

[5] Greydanus, S. et al. (2019). Hamiltonian Neural Networks. NeurIPS 2019.

[6] Hairer, E. et al. (2003). Geometric Numerical Integration. Springer.

[7] Kingma, D.P. & Ba, J. (2015). Adam: A Method for Stochastic Optimization. ICLR.